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ABSTRACT 


This  report  presents  new  and  recently  developed  concepts  which  are  useful 
for  obtaining  and  solving  equations  of  motion  of  multibody  mechanical  systems 
with  translation  between  the  respective  bodies  of  the  system.  These  concepts 
are  then  applied  in  the  study  of  human  head/neck  systems  in  high  acceleration 
conf lgurations . 

The  developed  concepts  include  the  use  of  Euler  parameters,  Lagrange's 
form  of  d'Alembert's  principle,  quasi-coordinates,  relative  coordinates,  and 
body  connection  arrays.  This  leads  to  the  development  of  efficient  computer 
algorithms  for  the  coefficients  of  the  equations  of  motion.  The  developed 
procedures  are  applicable  to  "chain-link"  systems  such  as  finite-segment  cable 
models,  mechanisms,  manipulators,  robots,  and  human  body  models. 

The  application  with  human  head/neck  models  consists  of  a 54  degree  of 
freedom,  three-dimensional  system  representing  the  head,  the  vertebrae,  and 
the  connecting  discs,  muscles,  and  ligaments.  The  computer  results  for  the 
system  in  a high  acceleration  configuration  agree  very  closely  with  available 
experimental  data. 
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INTRODUCTION  ! 

i 

Recently  there  has  been  considerable  Interest  in  the  development  of 
equations  of  motion  for  multi-body  mechanical  systems — that  is,  systems 
containing  many  rigid  bodies.  There  are  two  reasons  for  this  interest: 

First,  many  complicated  mechanical  systems  and  devices  such  as  manipulators, 
robots,  and  biosystems,  can  be  effectively  modelled  by  systems  of  rigid  bodies; 
and  second,  it  has  just  recently  been  possible,  with  the  aid  of  high-speed 
digital  computers,  to  obtain  efficient  numerical  solutions  of  the  governing 
dynamical  equations.  The  emphasis  of  researchers  working  with  multi-body 
systems  has  therefore  been  the  formulation  of  equations  of  motion  which  can 
easily  be  developed  into  numerical  algorithms  for  computer  codes. 

Most  of  this  recent  research  interest  has  been  with  multibody  systems 
consisting  of  linked  rigid  bodies  - that  is,  systems  of  connected  rigid  bodies 
such  that  adjacent  bodies  share  at  least  one  common  point  and  such  that  no 
closed  loops  or  circuits  are  formed.  Such  systems  are  sometimes  called 
"general-chain",  "open-chain",  or  "chain-link"  systems.  Figure  1.  depicts 
such  a system.  General  chain  systems  are  useful  for  modelling  chains,  cables, 
manipulators,  teleoperators,  antennas,  and  beams. 

There  are  some  systems,  however,  where  the  restriction  to  linked  rigid 
bodies  precludes  a satisfactory  modelling.  For  example,  with  a human  body 
model  it  is  frequently  advantageous  to  simulate  neck  stretch  during  periods 
of  high  acceleration  such  as  in  crash  environments.  Such  a simulation  is  not 
possible  with  a fully  linked  model.  Therefore,  it  is  of  interest  to  generalize 


the  multibody  models  to  Include  translation  between  the  bodies.  Figure  2. 
depicts  such  a generalization  of  the  system  of  Figure  1. 

This  report  presents  the  results  of  recent  research  efforts  to  develop 
efficient,  computer-oriented  algorithms  for  obtaining  and  solving  the  gover- 
ning dynamical  equations  of  motion  for  these  generalized  multibody  systems. 
The  report  also  contains  a summary  of  results  of  the  application  of  these 
procedures  with  human  head-neck  systems  in  high  acceleration  configurations. 

The  balance  of  the  report  is  divided  into  six  parts  with  the  first  part 
providing  a sunmary  of  earlier  efforts  to  model  multibody  systems.  This  is 
followed  by  two  parts  which  contain  the  general  geometrical  and  kinematical 
background  necessary  for  the  development  of  the  governing  equations.  The 
governing  equations  themselves  are  developed  in  the  next  part,  and  an  applica- 
tion of  the  developed  procedures  in  studying  head-neck  dynamics  is  presented 
in  the  subsequent  part.  The  final  part  contains  a summary  discussion  and 
suggestions  for  other  applications  of  the  developed  procedures. 
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PREVIOUS  MULTIBODY  SIMULATION  EFFORTS 


References  [1-36]*  provide  a summary  of  approaches  taken  to  obtain 
efficient,  computer-oriented  formulation  of  the  equations  of  motion  of 
multibody  systems  such  as  In  Figure  1.  In  one  of  these  approaches  [19,29,33], 
It  is  shown  that  it  is  possible  to  obtain  expressions  for  the  governing 
equations  In  a form  where  the  coefficients  are  easily  evaluated  through 
computer  algorithms.  This  approach  uses  Lagrange's  form  of  d'Alembert's 
principle,  as  exposited  by  Kane  and  others  [37,38,39],  together  with  relative 
orientation  coordinates  [40,41,42],  to  obtain  the  governing  equations. 
Although  this  principle  is  not  as  widely  used  as,  for  example,  Newton's 
laws  or  Lagrange's  equations,  it  has  the  advantage  of  automatic  elimination 
of  non-working  internal  constraint  forces  without  the  introduction  of  tedious 
differentiation  or  other  calculations. 

Recently,  it  has  been  suggested  by  Huston,  et.al.,  [42,43],  that  further 
efficiencies  in  the  development  and  solution  of  the  governing  equations  could 
be  obtained  through  the  use  of  Euler  parameters  as  described  by  Wittaker  [44] 
and  Kane  and  Likins  [45],  together  with  the  quasi-coordinates  suggested  by 
Kane  and  Wang  [46].  Specifically,  it  is  claimed  [42,43]  that  using  Euler 
parameters  together  with  relative  angular  velocity  components  as  generalized 
coordinate  derivatives  allows  for  the  avoidance  of  geometrical  singularities 
encountered  with  using  Euler  angles  or  dextral  orientation  angles  to  define 
the  relative  orientation  of  bodies.  (Recall  that  Euler  angles  may  be  defined 
by  aligning  mutually  perpendicular  axes  fixed  in  the  bodies  and  then 
successively  rotating  one  body  relative  to  the  other  about  the  third,  first. 
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; and  third  axes,  whereas  dextral  orientation  angles  may  be  defined  by 

successive  rotations  about  the  first,  second,  and  third  axes.) 
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PRELIMINARY  GEOMETRICAL  CONSIDERATIONS 

Body  Connection  Array 

Consider  a mechanical  system  such  as  depicted  in  Figure  1.  To 
develop  an  accounting  routine  for  the  system's  geometry,  arbitrarily 
select  one  of  the  bodies  as  a reference  body  and  call  it  B^.  Next, 
number  or  label  the  other  bodies  of  the  system  in  ascending  progression 
away  from  as  shown  in  Figure  1.  Now,  although  this  numbering  procedure 
does  not  lead  to  a unique  labeling  of  the  bodies,  it  can  nevertheless  be 
used  to  describe  the  chain  structure  or  topology  through  the  "body  connection 
array"  as  follows:  Let  L(k),  k*l,...,N  be  an  array  of  the  adjoining 
lower  numbered  body  of  body  B^,  For  example,  for  the  system  shown  in 
Figure  1.,  L(k)  is: 

L(k)  - (0, 1,1,3, 1,5, 6, 7, 6)  Cl) 

where 

(k)  - (1,2, 3, 4, 5, 6, 7, 8, 9)  (2) 

and  where  0 refers  to  an  inertial  reference  frame  R.  It  is  not  difficult 
to  see  that,  given  L(k),  one  could  readily  describe  the  topology  of  the 
system.  That  is.  Figure  1.  could  be  drawn  by  simply  knowing  L(k).  It  is 
shown  in  the  sequel  that  L(k)  is  useful  in  the  development  of  expressions  of 
klnematlcal  quantities  needed  for  analysis  of  the  system's  dynamics. 
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Transformation  Matrices 

Next,  consider  a typical  pair  of  adjoining  bodies  such,  as  B.  and  B 

J 21 

as  shown  in  Figure  3.  The  general  orientation  of  B^  relative  to  B^ 
may  be  defined  in  terms  of  the  relative  orientation  of  the  dextral 
orthogonal  unit  vector  sets  and  n^  Ci“l»2,3)  fixed  in  B^  and  B^ 
as  shown  in  Figure  2.  Specifically  n^  and  n^f  are  related  to  each 
other  as 

Sji  ■ Oi 

where  SJK  is  a 3 x 3 orthogonal  transformation  matrix  defined  as  [47J: 


* 5ji  • Sta 


CRegardlng  notation,  the  J and  K in  SJK.  and  the  first  subscripts  on  the 
unit  vectors  refer  to  bodies  B^  and  B^,  and  repeated  indices,  such  as  the 
m.  In  Equation  C3)  signify  a sum  over  the  range  Ceg.  1,...,3)  of  that 
Index.  Thus,  with  a computer  SJK^  would  be  the  array  SJK(J,M).) 

From  Equation  (.3) , it  is  easily  seen  that  with  three  bodies  B^ , Bfc, 

B. , the  transformation  matrix  obeys  the  following  chain  and  identity  rules: 


SJL  - SJK  SKL 
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and 


SJJ  - I - SJK  SKJ  - SJK  SJK 


-1 


(6) 


where  I is  the  identity  matrix. 

These  expressions  allow  for  the  transformation  of  components  of 
vectors  referred  to  one  body  of  the  system  into  components  referred 
to  any  other  body  of  the  system  and,  in  particular,  to  the  inertial  . 
reference  frame,  R.  For  example,  if  a typical  vector,  V,  is  expressed 
as 


_ (k)  „ (0) 

y*v  \i* Y n 


“0i 


(7) 


then 


V1(0)  - SOKjj  Vj 


(8) 


where  0 refers  to  the  inertial  frame,  R. 

Since  these  transformation  matrices  play  a central  role  throughout 
the  analysis,  it  is  helpful  to  also  have  an  algorithm  for  their  derivative. 


especially  the  derivative  of  SOK.  Using  Equation  (3) , and  noting  that  n^ 


are  fixed  in  R,  the  following  is  obtained: 


dCSOK^/dt  - nQ1  • “d 


/dt 


(9) 
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where  the  R in  d n^/dt  indicates  that  the  derivative  is  computed  in  R. 
However,  since  the  n^  are  fixed  in  B^,  their  derivatives  may  be  written 
as  o>k  x n^  where  is  the  angular  velocity  of  in  R.  Equation  (9) 


may  then  be  written  as: 


d(SOKy)/dt  - ^ 


or  as 


d(SOK) /dt  - WOK  SOK 


where  WOK  is  a matrix  defined  as 


WOK.  - -e 


im  eimn<likn 


and  where  to.  are  the  components  of  u.  referred  to  n and  e.  is  the 
kn  _k  .on  imn 

standard  permutation  symbol  [47,48],  (WOK  is  simply  the  matrix  whose  dual 
vector  [48]  is  <o . ) Equation  (11)  thus  shows  that  the  transformation  matrix 
derivative  may  be  computed  by  a simple  matrix  multiplication. 


Euler  Parameters 


Finally,  consider  describing  the  relative  orientation  of  Bj  and  B^ 
by  using  the  so-called  Euler  parameters  as  discussed  by  Whittaker  [44]  and 
Kane  and  Likins  [45].  It  is  well  known  [44]  that  B^  may  be  brought  into  any 
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general  orientation  relative  to  by  means  of  a single  rotation  about  an 
appropriate  axis.  If  X^  is  a unit  vector  along  this  axis  and  if  9^  is  the 
rotation  angle,  the  four  Euler  parameters  describing  the  orientation  of  B^ 
relative  to  B^  may  be  defined  as: 


eki  ‘ \i  >ln<v2) 


ek2  ‘ Xk2  sln(V2) 


ek3  ‘ lk3  5to<V2) 


Sa  - cos  Cek/2) 


where  the  X^  (i-1,2,3)  are  the  components  of  Xfc  referred  to  n^,  the  unit 
vector  fixed  in  B^.  Clearly,  the  (i**l,2,3,4)  are  not  independent  since: 


These  parameters  may  be  related  to  angular  velocity  components  by 
using  the  transformation  matricies  as  follows:  It  is  shown  in  [44,45]  that 
SJK  may  be  expressed  in  terms  of  these  parameters  as: 
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SJK 


£kl"£k2_£k3+£k4  2(eUek2"ek3sk4>  2<£kl£k3+£k2£k4) 


2<*kl£k2+£k3*k4> 


2 2 2 . 2 
“ckl+£k2'£k3+£k4 


2(ek2ek3-'kl£k4> 


2(8kl£k3-£k2sk4> 


2(£k2£k3+£kl£k4> 


2 2,22 

"ekl"ek2+ek3  ek4 


(15) 


Now,  by  solving  Equations  (11)  and  (12)  for  the  angular  velocity  components, 
one  obtains: 


<0^  - SOi^  S0K31  + SOK^  SOK32  + SOK23  SOK33 


“k2  - S0K31  S0K11  + SOK32  S0K12  + SOK33  S0K13 


(16) 


0^3  - S0K11  S0K21  + S0K12  S0K22  + SOK^  SOK23 


where  the  dot  designates  time  differentiation.  By  using  Equation  (15) , 


these  expressions  may  be  used  to  express  the  n. . components  of  the  angular 


velocity  of  B^  relative  to  B^  in  terms  of  the  Euler  parameters  as: 


\l  " 2(£k4  Skl  " £k3  S2  + £k2  S3  “ £kl  £k4) 


“k2  " 2<£k3  Si  + S4  S2  " £kl  ®k3  - S2  S4> 


(17) 


“k3  " 2(“£k2  Si  + Si  S2  + ek4  S3  " S3  S45 


.L.4H-  - X-7L*»^  ■ 


sq 
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(Regarding  notaCion,  in  the  sequel  "hats"  refer  to  relative  angular 
velocity  vectors  or  their  components.  That  is, the  represent  the 
angular  velocity  of  in  R and  uj^  represent  the  angular  velocity  of 
B^  relative  to  B^,  its  adjoining  lower  numbered  body.)  Equation  (17) 

A ^ 

may  now  be  solved  for  the  (i**l,...,4)  in  terms  of  the  leading 

to  the  expressions: 


£kl  “ ^S 4 “kl  + ek3  “k2  " ek2  “k3) 

Ek2  " «“S3  Si  + Sk4  “k 2 + Si  “k3} 

(18) 

ek3  “ Jj(ek2  Si  " £kl  “k2  + Sk4  “k3) 
ek4  ” ^"Sl  Si  " ek2  “k2  " S3  S3} 


This  solution  is  quickly  obtained  by  observing  that  if  Equation  (14)  is 
differentiated  and  placed  with  Equation  (17),  the  resulting  set  of  equations 
could  be  written  in  the  matrix  form: 


* 

Si 

S4 

"ek3 

Ek2 

"Si 

A 

“k2 

- 2 

S3 

Ek4 

"Si 

-S2 

S3 

-Ck2 

£kl 

S4 

"S3 

A 

\4 

Si 

Sk2 

S3 

S4 

— — 

(19) 
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where  Is  equal  to  the  derivative  of  Equation  (14)  and  has  the  value 
zero.  The  square  matrix  In  Equation  (19)  is  seen  to  be  orthogonal 
(le.  the  Inverse  is  the  transpose)  and  hence.  Equations  (18)  follow 
immediately  from  (19)  upon  letting  be  zero. 
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KINEMATICS 

Coordinates 

A multibody  system  of  N bodies,  with  translation  permitted  between 
the  bodies  will,  in  general,  have  6N  degrees  of  freedom.  Let  these  be 
described  by  6N  generalized  coordinates  x^  (£**1, . . . ,6N)  and  let  the  first 
3N  of  these  be  divided  into  N triplets  describing  the  relative 
orientation  of  the  successive  bodies  of  the  system.  Let  the 
remaining  3N  x^  also  be  divided  into  N triplets  representing  the  relative 
displacement  of  the  successive  bodies  of  the  system.  As  before,  let 
be  a typical  body  of  the  system  and  let  be  its  adjacent  lower  numbered 
body,  as  in  Figure  3.  The  angular  velocity  of  B^  relative  to  B^  (that  is, 
the  relative  rate  of  change  of  orientation)  may  then  be  written  as: 

k • “u  5jl  + “k2  5J2  + “k3  5j3  <20) 

where  n^  (j*l,...,Nj  i*l,2,3)  are  mutually  perpendicular  dextral  unit 
vectors  fixed  in  B^.  Next,  let  these  bodies  be  displaced  relative  to 
each  other  with  the  displacement  measured  by  the  vector  as  shown  in 

Figure  4.,  where  0^  and  0^  are  arbitrarily  selected  reference  points  of 

Bj  and  B^.  0^,  which  is  fixed  in  B ^ , is  the  connection  point  or  "origin" 

of  B.  . Then  £.  may  be  written  in  the  form: 

JC 
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In  g«n«r<l,  Equations  (23)  are  non-integrable.  That  is,  they  cannot 
be  Integrated  to  obtain  generalized  orientation  coordinates  » 

x3k-1  , . Thus,  explicit  parameters  *3^2’  x3k-l*  and  x3k  do  not 

in  general  exist — hence,  the  name  "quasi-coordinates".  However,  since 
parameters  are  needed  to  relate  the  relative  orientation  of  the  bodies 
to  the  respective  relative  angular  velocities,  let  the  Euler  parameters 
introduced  in  the  foregoing  section  be  used  for  this  purpose.  Hence,  if  the 
orientation  of  a typical  body  relative  to  is  described  by  the  four 
parameters  (i-1,. . . ,4) , the  geometry  and  kinematics  of  the  entire  system 
may  be  expressed  in  terms  of  the  4N  Euler  parameters  (k»l,...,N;  i-1,..., 4), 
the  3N  relative  angular  velocity  components  (k-l,...,N;  i-1, 2, 3),  and  the 
3N  displacement  components  (k-1,... ,N;  i-1,2,3). 

Angular  Velocity 


I 


I 


1 

t 

i 
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The  angular  velocity  of  a typical  body  B^  in  the  inertial  frame  R is 
readily  obtained  by  the  addition  formula  as  [38): 

A ^ 

+ (25) 


where  the  relative  angular  velocities  on  the  right  side  of  this  expression 
are  each  with  respect  to  the  respective  adjacent  lower  numbered  bodies  and 
where  the  sum  is  taken  over  the  bodies  of  the  chain  from  outward  through 
the  branch  containing  B^.  The  L(k)  array  introduced  in  the  foregoing  section 
can  be  useful  in  computing  this  sum:  Consider  for  example,  the  system  shown 
in  Figure  1.  The  angular  velocity  of  B^  is: 
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“9  " “l  + "5  + “6  + “9 


The  subscript  Indices  (ie.  9, 6, 5,1)  may  be  obtained  from  L(k)  as 

follows:  Consider  L(k)  as  a function  mapping  the  (k)  array  (See  Equation  (2)) 

into  the  L(k)  array.  Then,  using  the  notation  that  L^(k)  - (k), 

LX(k)  - L(L°(k) , L2(k)  - L(LX(k),  LJ  (k)  - L(LJ_1(k)),  it  is  seen 

(see  Equation  (1))  that: 


L°(9)  - 9,  LX(9)  - 6,  L2(9)  - 5,  L3(9)  - 1 


Therefore,  u.  may  be  written  aa: 

•*7 


u9  - l w , q - LP(9) 


Hence,  in  general,  the  angular  velocity  of  may  be  written  as: 


ov  " l <d  , q - LP(k)  C29) 

^ p-0  ~q 

where  r is  the  index  such  that  Lr(k)  ■ 1 and  it  is  obtained  by  comparing 
LP(k)  to  1.  The  index  r represents  the  number  of  bodies  from  to  B^  in 

that  branch  of  the  chain  system  B^.  For  example,  for  the  system  of  Figure  1. , 
if  k*9,  r»3.  Equation  (29)  is  thus  an  algorithm  for  determining  once 
and  L(k)  are  known. 


16 


Ik  “ 5kl  Sjl  + *k2  2j2  + *k3  2j3 


(21) 


Following  Kane  and  Wang  [46] , Introduce  6N  parameters  (Jt*l, . . . ,6N) 
defined  as: 


y4  - x£  £ - 1,...,  6N 


where  the  first  3N  of  these  are 


(22) 
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By  examining  Equations  (20,  (23),  and  (25)  it  Is  seen  that  may  be 


written  in  the  form 


y.  n (30) 

-k  klm  - om 

where  there  is  a sum  over  the  repeated  indices  and  where  (k-l,...,N; 

fc*l,...,3H;  m-1,2,3)  form  a block  array  of  coefficients  needed  to  express 

to,  in  terms  of  n . In  view  of  Equations  (3),  (16),  (20),  and  (23), 

~k  -om 

it  is  seen  that  the  elements  of  the  array  may  be  obtained  from  the  SOK 

transformation  matrices.  Moreover,  it  can  be  shown  that  the  matching  between 
the  elements  of  the  w.  ^ and  SOK  arrays  is  solely  dependent  upon  the  body 
connection  array  L(k). 


To  see  this,  consider  for  example  the  angular  velocity  of  of  the 
system  of  Figure  1:  From  Equation  (25),  is 


AAA 


0)  « a»  + u>_  + 0). 
~4  -**4 


where  from  Equations  (3),  (20),  and  (23)  5^,  u^,  and  may  be  written 


Sl  ■ yl2oi  + y2  502  + y3  2(33  ■ ri  Smj  2c 


S3  ■ y7  Su  + yS  Sl2  + y9  513  ‘ ywj  S01.h  2c 


24  ■ y10  231  + yll  232  + y12  233  * Vj  S03mJ  2o» 


«»£  1.2.3 


0 i - 4,5,6 


u>4Am  - SOl^g  l - 7,8,9  m - 1,2,3  (35) 


S03mJl-9  1 * 10.H.12 


£ > 12 


where  6^  are  the  identity  matrix  components  [47,48]. 

Next,  consider  that  the  results  such  as  Equation  (35)  may  be  obtained  for 

the  entire  system  of  Figure  1.  or  Figure  2.  from  a table  such  as  Table  1., 

where  the  "m"  entries  of  the  array  are  the  column  of  the  transformation 

matrices.  Finally,  note  that  the  non-zero  entries  in  a typical  row, 

say  the  kth  row  of  Table  1.  are  obtained  as  follows:  Let  P * L(k). 

til  * 

Then  SOP  is  placed  in  the  k column  of  triplets  of  x^.  Next,  let  Q»L(P). 

til  • 

The  SOQ  is  placed  in  the  P column  to  triplets  of  x^,  etc.  That  is,  SOM 
is  placed  in  column  L^  ^(k)  where  M ■ L^(k),  j*l,...,r+l  with  r determined 
from  Lr(k)  ■ 1. 

Finally,  it  is  interesting  to  note  that  the  elements  of  the  10^^  array 
(and  hence,  the  transformation  matrix  columns  of  Table  1.)  are  components 
of  the  "partial  rate  of  change  of  orientation  vectors"  as  originally  defined 
by  Kane  [37] . 


* 


I 

I 

I 

I 

I 

I 
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Angular  Acceleration 

The  angular  acceleration  of  in  R may  be  obtained  by  differentiating 
Equation  (30).  Noting  that  the  n are  constant,  this  leads  to: 


2k  " (“kl»  y + V.  V 2am 


(36) 


f 

! 

* 


5 


* 


r 


j 


A table  containing  the  can  be  constructed  directly  form  the  corresponding 

table  for  the  For  example,  for  the  system  of  Figure  1.,  such  a 

table  is  shown  in  Table  2. 

Mass  Center  Velocities 

The  velocity  and  acceleration  of  the  mass  center  of  a typical  body 
(k-l,...,N)  may  be  obtained  as  follows:  Let  r^  locate  G^  relative 
as  shown  in  Figure  4.  Since  0^  is  located  relative  to  by  and  if  is 
located  relative  to  0.  by  the  vector  q,  (See  Figure  4.),  then  by  continuing 

j 

this  procedure,  G^  may  ultimately  be  located  relative  to  a fixed  point  0 in 
R,  the  inertial  reference  frame.  For  example,  for  Body  Bg  of  Figure  2.,  the 
position  vector  Pg  of  Gg  relative  to  0 is: 

Is  " k + 35  + Is  + 36  + 5e  + 3?  + I7  + 3s  + Is  + Is  (37) 

In  general,  for  Body  B^,  the  position  vector  of  relative  to  0 is: 


A 


4 


A 
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where  s ■ L^(k),  S ” L**+^(k),  and  u is  the  index  such  that  Lu(k)  ■ 1, 
and  where  q^  is  0.  By  differentiating,  the  velocity  of  in  R is  obtained 


Zk  - (s6Kih  rkh  + \ ts6sih  <’.h  + W 
q*0 
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where  WK  . is  defined  as: 
mh 


K . 9 - =*  SOK  - 

327z  ph 


smpi  “kji  S0Kph 


v^(3jj+-£)m  * ®<“1»»»»»3N;  ®*1»2,3) 


Mass  Center  Accelerations 


Similarly,  by  differentiation  of  Equations  (40),  the  acceleration 
of  in  R is 


Sk  ' <Vm  yl  + vUo  yl>5c 


where  the  non-zero  are,  by  Equations  (41)  to  (43), 


Vk£m  " ^mhJl  rkh  + ^ ^W^mh2  ^sh  0^1»***»N;  i-l,...,3N,  m-1,2,3) 

q-0 


+ qsh  + WSmhS,  hh] 


where  is: 


"W  - -mpl  <Vl  S0Kph  + Vl  S°V 
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EQUATIONS  OF  MOTION 


Consider  again  a general  chain  system  such  as  shown  in  Figure  2.,  and 
imagine  the  system  to  be  subjected  to  an  externally  applied  force  field. 

Let  the  force  field  on  a typical  body  B^,  be  replaced  by  an  equivalent 
force  field  consisting  of  a single  force  F.  , passing  through  G.  together 

*»lC  tC 

with  a couple  with  torque  M^.  Then  Lagrange's  form  of  d'Alembert's  principle 
leads  to  governing  dynamical  equations  of  motion  of  the  form  [38]: 


+ Fl*  " 0 1 “ 1 6N 


F^  (2,-1,... ,6N)  is  called  the  generalized  active  force  and  is  given 


FA  " Vk2m  Fkm  + <\£m  **km 


where  there  is  a sum  from  1 to  N on  k and  from  1 to  3 on  m,  and  where  F 

* n 

and  are  the  components  of  F^  and  with  respect  to  nQn  . F^* 

(2.-1,. . . ,6N)  is  called  the  generalized  inertia  force  and  is  given  by: 


V - vk*m  Fk£  + U “fcn 
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where  Che  Indices  follow  the  same  rules  as  in  Equation  (48),  and  where 
F^£  and  are  n^  components  of  inertia  forces,  F^*,  and  inertia  torques, 
M^*,  given  by  [38] : 


Ffc*  - -n^  ^ (no  sum) 


“k*  * "in  * 5k  “ “k  x (h  * ^ (no  sua)  (52) 

where  m^  is  the  mass  of  and  1^  is  the  inertia  dyadic  of 
relative  to  Gfc  (k-l,...,N).  (F*,  with  line  of  action  passing  through 

together  with  M*  are  equivalent  to  tlie  inertia  forces  on  B^  [38].) 
Through  use  of  the  shifter  transformation  matrices,  1^  may  be  written 
in  the  form: 

I.  - L n n (53) 

~k  Kmn  ~om  -.on  7 


By  substituting  Equations  (36)  and  (44)  into  Equations  (51)  and  (52) 
and  ultimately  into  Equation  (47) , the  equations  of  motion  may  be  written 
in  the  form: 
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aJlp  yp  “ *1 


where  there  is  a sum  from  1 to  6N  on  p and  where  and  f ^ are 
given  by: 


lJLp  ^ vkpm  vkAm  + *kmn  ^kpm 


m + “k  Vkim  vkqm  yq  + ^kmn  ^Hcim  ^kqn  yq 


+ enmh  *kmr  ^kqn  ^ksr  ^kin  yq  ys^ 


where  there  is  a sum  from  1 to  N on  k,  from  1 to  6N  on  q and  s,  and  from 

1 to  3 on  the  other  repeated  indices. 

Recall  that  the  first  3N  y are  relative  angular  velocity  components. 

P 

These  may  be  related  to  the  Euler  parameters  by  N sets  of  first  order 

equations  of  the  form  of  Equations  (18). 

Equations  (54) , (20) , and  the  4N  equations  of  the  form  of  Equations 

(18)  form  a set  of  13N  simultaneous  first-order  differential  equations  for 

the  6N  y^,  the  3N  and  the  4N  Euler  parameters  (h«l,...,N; 

1*1,... ,4).  Since  the  coefficients  a^  and  f^  in  Equations  (54)  are 

algebraic  functions  of  the  physical  parameters  and  the  four  block  arrays 
• • 

^klm'  UkAm’  Vklm  an<*  Vk!m'  comPuter  algorithms  can  be  written  for  the 
numerical  development  of  these  governing  equations.  Moreover,  once  these 


* 
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arrays  are  developed,  the  system  of  equations  consisting  of  Equations  (54) , 
(20) , and  4N  equations  of  the  form  of  Equations  (18) , may  also  be  solved 
numerically  by  using  one  of  the  standard  numerical  integration  routines 
and  a linear  equation  solver. 

The  development  of  these  computer  algorithms  and  the  numerical  development 
of  Equations  (54)  might  proceed  as  follows:  First,  let  the  body  connection 
array  L(k)  (See  Equation  (1))  together  with  the  geometrical  and  physical 


parameters  r^»  g^,  I^»  and  m^  (See  Equations  (38),  (51),  and  (52).)  and 
the  applied  forces  and  moments  and  (See  Equation  (48).)  be  read  into 


the  computer.  (Let  r^,  g^,  1^  and,  if  desired,  F^  and  be  expressed  in 


terms  of  n^ . ) Next,  from  assumed  initial  values  of  form  the 


transformation  matrix  arrays  SOK  using  Equations  (15)  and  (5).  Use  these 


arrays  to  express  rfc,  gfc,  and  possibly  F^  and  in  terms  of  n^.  Next, 
using  L(k)  and  SOK  writs  an  algorithm,  with  Tables  1.  and  2.  as  a guide,  to 
form  and  • For  example,  to  obtain  the  non-zero  observe  that 

if  L(k)  “ p,  then  ■ SOP^  (m-1,2,3;  i-3p+l,  3p+2,  3p+3).  Then,  if 
L(p)  - q,  L2(k)  - q and  <\9m  ■ SOQ^  (m-1,2,3;  2-3q+l,  3q+2,  3q+3). 

This  assignment  procedure  is  continued  until  unity  is  reached  or  r times 
where  r is  given  by  Lr(k)  - 1 (See  the  remark  following  Equation  (29).). 
v^9m  and  nay  then  be  obtained  using  Equations  (40)  to  (47).  Finally, 
numerical  values  of  the  coefficients  a^  and  f^  of  the  governing  differential 
equations  (54)  may  then  be  obtained  from  Equations  (55)  and  (56).  These 


equations  may  then  be  integrated  numerically  to  obtain  incremental  values  to  the 

initial  values  of  the  parameters  y , £,, , and  x (p-1, . . . ,3N+3;  k-l,...,N;  1-1, 2, 3, 4 

p id  q 

and  q-1,2,3),  at  the  end  of  a time  interval,  say  t, . New  values  of  the 


repeated  until  a history  of  the  configuration  and  motion  of  the  system  is 
determined. 


Specific  computer  algorithms  following  this  general  procedure  have 
been  written  and  validated.  A listing  together  with  a tape  copy  (or  card 
deck)  are  available  at  reproduction  cost  from  the  authors. 
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APPLICATION  WITH  HEAD-NECK  SYSTEMS 

Previous  Simulation  Efforts 

Recently,  there  has  been  considerable  interest  in  using  the  foregoing 
procedures  in  the  modelling  of  biodynamic  systems.  Specifically,  there  has 
been  Interest  in  modelling  the  human  body  - and  particularly,  head-neck  systems  - 
during  periods  of  high  acceleration,  as  experienced  in  vehicle  accidents.  This 
interest  stems  from  the  fact  that  accident  injuries,  including  both  direct  and 
Indirect  (for  example,  "whiplash")  impact,  are  basically  mechanical  phenomena. 

The  emphasis  on  modelling  the  head-neck  system  is  stimulated  by  the  belief 
that  as  many  as  60  - 70%  of  vehicle  related  accident  fatalities  are  a direct 
result  of  injuries  to  the  head-neck  system. 

There  are  a number  of  head-neck  simulation  models  discussed  in  the  tech- 
nical literature.  Specifically,  in  1971,  Orne  and  Liu  [60]  developed  a discrete- 
parameter  spine  model  which  simultaneously  accounts  for  axial,  shear,  and 
bending  deformation  of  the  discs,  for  the  variable  size  and  mass  of  the  vertebrae 
and  discs,  and  for  the  natural  curvature  of  the  spine.  They  also  present  an 
extensive  literature  review  of  spine  models  prior  to  1970.  Later  in  1971, 
McKenzie  and  Williams  [61]  used  the  Ome-Liu  model  to  develop  a two-dimensional 
discrete-parameter  head -neck-tor so  model  for  "whiplash"  investigation.  A two- 
dimensional  mechanical  linkage  model  simulating  head-neck  response  to  frontal 
impact  has  been  presented  by  Becker  [62].  This  model  allows  for  elongation  of 
the  neck.  It  concentrates  the  mass  at  the  head  mass  center.  Springs  and 
dampers  are  used  to  control  the  elongation  of  the  model.  A three-dimensional 
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neck-torso  linkage  vehicle-occupant  model  has  been  developed  by  Bowman  and 
Robbins  [63].  The  model  has  two  ball-and-socket  joints  and  the  neck  can 
elongate  with  the  motion  limited  by  joint  stopping  moments. 


In  addition  to  these  computer  models,  there  have  also  been  developed  a 
number  of  anthropometric  dummy  models.  (These  are  currently  used  extensively 
by  the  automotive  industry.)  In  1972,  Melvin,  et.al.  [64]  presented  a mech- 
anical neck  for  authropometric  dummies . The  neck  consists  of  three  steel  uni- 
versal joints  pinned  into  aluminum  discs  with  shaped  rubber  discs  around  the 
joints.  The  joints  allow  the  neck  to  move  in  flexion,  extension,  and  lateral 
flexion  but  do  not  allow  for  either  rotation  or  elongation.  A mechanical  neck 
has  also  been  presented  by  Culver,  et.al.  [65].  It  consists  of  four  ball-joint 
segments  and  one  pin-connected  "nodding”  segment.  Viscoelastic  resistive 
elements  inserted  between  the  segments  provide  for  bending  resistance  and 
energy  dissipation  with  the  primary  objective  being  to  model  flexion  and  ex- 
tension responses. 

In  this  part  of  the  report,  there  is  presented,  as  an  application  of  the 
foregoing  procedures,  a comprehensive,  three-dimensional,  head-neck  computer 
model  which  has  54  degrees  of  freedom  and  includes  the  effects  of  discs, 
muscles,  and  ligaments.  The  model  is  developed  by  considering  the  skull  and 
vertebrae  as  a chain  system  of  rigid  bodies  which  may  translate  relative  to 
one  another.  The  soft  tissue  effects  of  the  discs,  muscles,  and  ligaments 
are  modelled  by  nonlinear  springs  and  dampers  between  the  bodies.  The  model 
is  based  primarily  on  the  research  of  J.  Huston  and  Advani  [55,56,57], 


The  balance  of  this  part  of  the  report  contains  a description  of  the 
modelling  itself  and  the  development  of  the  governing  dynamical  equations  of 
motion.  This  is  followed  by  a comparison  of  results  from  numerical  integration 
of  these  equations,  with  available  experimental  data. 


Head-Neck  Mod el line 


A comprehensive  presentation  of  the  head-neck  anatomy  may  be  found  in 
references  [66-73],  The  anatomy  is  conveniently  divided  into  two  categories: 
bones  and  soft  tissue. 


Bones 


The  largest  and  heaviest  is  the  skull  which  consists  of  a large  cranial 
cavity  (enclosing  the  brain)  and  smaller  bones  of  the  face  and  jaw.  The  skull 
is  actually  composed  of  21  closely  fitted  bones.  The  other  bones  of  the  head- 
neck  system  are  seven  cervical  vertebrae  (C1-C7)  which  support  and  provide 


mobility  to  the  head.  The  first  of  these  Cl,  called  the  "atlas",  supports  the 
skull.  The  second  C2,  called  the  "axis",  is  distinctive  because  of  its  adontoid 
process  (or  axis)  which  rises  perpendicularly  to  the  vertebrae.  The  five 
remaining  cervical  vertebrae  are  roughly  annular  in  shape  and  are  similar  to 
each  other  with  a slight  increase  in  size  going  down  from  C3  to  C7. 


Soft  Tissue 


The  soft  tissue  is  composed  primarily  of  the  discs,  the  muscles,  the 
ligaments,  and  the  brain.  The  discs  provide  the  cushioning  or  separation  for 
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Che  vertebrae.  They  are  annular  In  shape.  The  ligaments  connect  the  cervical 
vertebrae  to  each  other  and  thus  allow  for  the  gross  and  fine  movement  of  the 
head  and  neck.  The  muscles  control  the  movement  of  the  head  and  neck  which  may 
be  classified  grossly  as:  flexion,  extension,  and  rotation.  The  muscles 
originate  on  the  various  cervical  vertebrae,  the  skull,  the  spine,  and  the 
shoulder  bones.  The  brain  tissue  is  basically  four  mass  volumes  composed  of 
two  cerebral  hemispheres  In  the  upper  half  of  the  skull,  the  triangular  shaped 
cerebellum  In  the  lower  posterior  and  the  brain  stem  in  the  center  of  the  skull. 


Modelling 


The  head-neck  system  is  modelled  by  a system  of  9 rigid  bodies  representing 
the  skull,  vertebrae,  and  torso  as  shown  in  Figure  5.  and  springs  and  dampers 
representing  the  discs,  ligaments,  and  muscles.  The  masses,  inertia  matrices, 
and  overall  geometry  of  the  rigid  bodies  are  adjusted  to  match  the  actual 
human  values  [70].  Each  body  has  6 degrees  of  freedom  and  hence,  the  entire 
system  has  a total  of  54  degrees  of  freedom. 

Following  Orne  and  Liu  [60]  the  discs  are  modelled  in  the  axial  direction 
as  two-parameter  viscoelastic  solids  with  the  uniaxial  force-displacement 


relationship  being: 


F - (A/h) (d16  + d26) 


In  bending  and  shear  the  discs  are  modelled  as  linear  elastic  solids.  Using 
the  principles  of  strength  of  materials  theory  [70],  the  following  force  and 
moment  equations  are  developed: 


J 


(59) 


where  P^  and  are: 


and 


Fy  - ( -r  > < -T  + 9x)/P2 


**-(£>  (dl  5z  + d2<Sz) 


M - 


El,  66 

< hPT  > I -h*  + <P2+3>  9x  > 


El,  -66 


M 


* hP2  ^ * 


— + (Pl+3)  9y  1 


M - JG0  /h 
z z 


P1  * 1 + 


12E1lkx 


GAh 


12EI.k 

1+ k* 


GAh 


(60) 


(61) 


(62) 


(63) 


(64) 


(65) 


where  as  shown  in  Figure  5.,  Z is  in  the  axial  (up  direction,  X is  forward  and 


Y is  to  the  left 
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The  ligaments  are  modelled  as  non-linear  elastic  bands  capable  of  exerting 
force  only  in  tension.  The  force-displacement  relation  is  taken  as: 


F - JLjd  + l2$ 


The  muscles  are  modelled  as  two-parameter,  visco-elastic  solids,  which, 
like  the  ligaments,  only  exert  force  when  in  tension.  The  force-displacement 
relation  is  taken  as: 


F ® m^6  + 


The  joint  constraints  (limiting  the  relative  motion  of  the  bodies)  are 
modelled  as  one-way  dampers.  The  force-displacement  and  moment-rotation  relation 
are  taken  as: 


-cS  for  5>0 


0 for  5<0 


— C§  for  8>0 


0 for  9<0 


where  the  damping  constant  is 


C + C,  (X-Xmax)  for  X > Xinax 
o 1 


C - C 


for  Xmin  < X < Xmax 


C + C, (X-Xmin)  for  X < Xmin 
o l 


where  X,Xmax,  and  Xmin  are  the  values  of  the  displacement  or  rotations  variable 
and  its  corresponding  maximum  and  minimum  values. 


The  values  of  these  various  constants  for  the  discs,  ligaments,  muscles, 
and  joints  for  the  various  directions  and  motion  are  difficult  to  specify  pre- 
cisely due  to  a lack  of  experimental  data.  However,  the  values  for  the  discs 

■ may  be  obtained  from  Markold  and  Steidal  [74],  Ome  and  Liu  [60],  and  McKenzie 

and  Williams  [6]].  The  ligament  and  muscle  attach  points  may  be  obtained  from 
Francis  [75],  Lanier  [76],  and  Todd  and  Lindala  [77],  with  the  spring  and  vis- 
coelastic constants  obtained  from  Nunley  [78]  and  Close  [79]. 

I 

| 

Governing  Equations 

l 

\ 

' 

Tht  procedures  developed  in  the  foregoing  parts  of  the  report  are  directly 
applicable  to  the  model  of  Figure  5.  including  the  simulated  disc,  muscle,  and 

■ ligament  forces.  Specifically,  as  noted  earlier,  the  model  has  54  degrees  of 
freedom  (27  translation  and  27  rotation) . This  leads  to  a system  of  117  simul- 
taneous first-order  differential  equations  of  the  form  of  Equations  (18) , (20) , 
and  (54).  The  disc,  muscle,  and  ligament  forces  are  included  in  the  generalized 
active  forces  F^  of  Equation  (56). 

Comparison  with  Experimental  Data 

It  is  difficult  to  obtain  experimental  data  which  is  suitable  for  checking 
the  model.  This  is  due  to  the  expense  and  impracticality  of  using  dummies. 


cadavers,  or  animal  surrogates  and  due  to  the  limited  experimental  range  with 
human  volunteers.  However,  several  experiments  have  been  conducted  which  may 
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be  used  to  obtain  a validation  of  the  model.  In  one  of  these,  a seated  cadaver 
was  subjected  to  head  impacts  by  a rigid  pendulum.  Accelerometers  were  used  tc 
measure  the  resultant  frontal  and  occipital  head  impact  forces  and  accelerations. 
Using  the  impact  force  data  as  input,  the  acceleration  was  calculated  using 
the  computer  model.  A comparison  of  the  results  for  two  of  the  frontal  impact 
experiments,  6-2  and  6-5  is  shown  in  Figures  6.-9. 

In  the  same  set  of  experiments,  high-speed  cameras  were  used  to  measure 
the  acceleration,  velocity,  and  displacement  of  the  mass  center.  A comparison 
of  the  results  with  those  predicted  by  the  computer  model  for  experiments  6-1 
and  6-2  are  shown  in  Figure  10. 

Finally,  the  model  was  checked  against  live  human  data  generated  by  Ewing 
and  Thomas  [33]  using  elaborate  testing  facilities.  A comparison  of  the  results 
for  the  head  angular  acceleration,  angular  velocity,  and  angular  displacement 
is  shown  in  Figures  11.,  12.,  and  13. 
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I DISCUSSION  AND  CONCLUSIONS 


The  results  of  using  the  modelling  procedures  outlined  herein  and  numerically 
integrating  the  resulting  governing  differential  equations  (54)  for  a number 
of  other  physical  systems  (in  addition  to  head-neck  systems)  are  reported  and 
discussed  in  References  [40,41,49,50,51,52,53]. 


The  application  of  Equations  (54)  with  these  systems,  however,  is  based 
on  the  use  of  relative  orientation  angles  between  the  respective  bodies  of  the 

i 

system  as  the  generalized  coordinates  (x^)  as  opposed  to  the  use  of  Euler 
parameters  and  quasi-coordinates  as  outlined  in  the  foregoing  sections.  A < 

problem  which  arises  in  the  numerical  solution  of  Equations  (54)  where  orienta- 
tion angles  are  used  is  that  there  always  exists  values  of  the  angles  and  hence,  • 

configurations  of  the  system,  for  which  the  determinant  of  a^  is  zero.  A 
numerical  solution  will,  of  course,  fail  to  converge  at  these  singular  configura- 
tions of  the  system,  and  convergence  is  very  slow  for  configurations  in  the 
vicinity  of  a singularity.  This  problem  is  avoided  by  using  Euler  parameters 
| to  relate  the  orientation  geometry  to  the  angular  velocity.  ] 


The  advantages  of  using  Lagrange's  form  of  d'Alembert's  principle  to 
obtain  the  governing  equations  of  motion  for  multi-body  mechanical  systems 
has  been  exposited  in  detail  in  References  [29]  and  [39].  Basically,  this 
principle  has  the  advantages  of  Lagrange's  equations  or  of  virtual  work  in 
that  non-working  internal  constraint  forces,  between  the  bodies  of  the  system, 
are  automatically  eliminated  from  the  analysis,  and  may  therefore  be  ignored 
in  the  formulation  of  the  governing  equations.  The  principle,  however,  has 
the  additional  advantage  of  avoiding  the  differentiation  of  scalar  energy 
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accelerations  are  performed  by  vector  cross  products  and  multiplication  algor- 
ithms — procedures  which  are  Ideally  suited  for  numerical  computation.  As 
with  Lagrange's  equations,  Lagrange's  form  of  d'Alembert's  principle  requires 
the  use  of  generalized  coordinates  to  define  the  system  geometry.  The  use  of 
Euler  parameters  to  avoid  problems  with  singularities,  as  discussed  above, 
leads  naturally  to  the  use  of  relative  angular  velocity  components  as  the 
generalized  coordinate  derivatives.  This  in  turn  leads  to  additional  compu- 
tational advantages  as  observed  by  Kane  and  Wang  [46]  and  Likins  [54].  Speci- 
fically, by  using  relative  angular  velocity  components  as  the  principle  para- 
meters of  the  analysis,  the  coefficient  matrices  in  the  governing  equations 
can  be  obtained  directly  from  the  body  connection  array  L(k)  (See  Tables  1. 
and  2.). 

The  use  of  "relative"  coordinates,  that  is,  angular  velocity  components 
of  the  bodies  with  respect  to  their  adjoining  bodies,  as  opposed  to  "absolute" 
coordinates,  that  is,  angular  velocity  components  in  inertial  space,  also  con- 
tributes to  the  computational  advantage.  In  applications  with  specific  geo- 
metrical configurations  [40,41,49-53],  it  is  seen  that  the  geometry  is  more 
easily  described  in  terms  of  relative  coordinates. 

Finally,  the  generalization  to  allow  translation  between  the  bodies  of 
the  system  makes  the  analysis  applicable  to  a much  broader  class  of  problems 
than  was  possible  with  those  previous  analyses  which  are  restricted  to  linked 
multibody  systems.  For  example,  with  the  head-neck  system,  the  use  of  trans- 
lation variables  between  the  vertebrae  is  necessary  to  obtain  a satisfactory 
model  of  the  system.  Moreover,  this  generalization  to  include  translation  is 
a natural  extension  of  the  analyses  of  [33,42,49,50,51]. 
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Regarding  the  application  to  the  head-neck  system.  Figures  7.  - 13.  show 
there  is  agreement  between  the  experimental  results  and  ♦‘hose  predicted  by  the 
computer  model.  This  is  indeed  encouraging  and  it  suggests  that  this  head-neck 
model  represents  one  of  the  most  sophisticated  models  available.  However,  more 
testing  and  refining  needs  to  be  done.  Specifically,  the  three-dimensional 
features  of  the  model  need  to  be  further  checked  with  experimental  data.  Also, 
better  experimental  values  for  the  soft  tissue  mechanical  properties  need  to 
be  obtained.  Finally,  the  effect  of  muscle  time  delay  needs  to  be  Incorporated 
into  the  model. 

Beyond  this,  as  Injury  criteria  becomes  better  established,  the  model  can 
serve  as  an  effective  and  economical  tool  for  predicting  injury  in  a variety 
of  high-acceleration/high-accident  configuration  environments.  It  could  then 
be  used  for  the  development  and  design  of  safety  and  restraining  devices. 

Finally,  the  entire  analysis  and  the  procedures  outlined  in  this  report 
are  developed  with  the  intent  of  obtaining  efficiencies  in  a computer  or  numeri- 
cally oriented  development  and  solution  of  the  governing  dynamical  equations  of 
large  multibody  systems.  As  such,  its  most  productive  application  is  likely 
to  be  with  systems  such  as  finite-segment  biodynamic  models,  chains,  cables, 

; robots,  manipulators,  teleoperators,  etc. 
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Figure  3.  Two  Typical  Adjoining  Bodies 
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Figure  10.  Comparison  of  Model  and  Experiment  for  Head  Mass-Center 
Displacement,  Velocity,  and  Acceleration 
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Table  2.  kim  for  the  System  of  Figure 
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